% *****************************************************************************
%%                                                                            *
%%                                                                            *
% AUTHOR: Vigneshwar Ramakrishnan											  *	
% Date created: 29 August 2022												  *	
% CC-BY-SA-NC                                                                 *
% Caculation of the angular velocity and rotational density of states		  *	
%%                                                                            * 
% COM and atomic velocity obtained from simulations. 						  *	
% The angular velocity is obtained as w = I^-1*L							  *	
% where L = sum(m_j*(r_vec_j x v_vec_j)) where j is the atom (H or O) and 	  *
% r_vec_j denotes the position vector of atom j w.r.t COM of molecule         *
% I(ii) = sum (mass*(r(i)^2), i = x,y,z										  *
%                                                                             *
%******************************************************************************

% Constants
clear;clc;
warning off all

kb = 1.3806503E-23;                 % boltzmann constant: m^2 kg s^-2 K^-1
T = 300;                           % Temperature: K
NA = 6.023E23;                     % Avogadro number
mwt_water = 0.018;                 % mol. weight of water: kg/mol (18 g/mol)
m_water = mwt_water/NA;            % mass of 1 water molecule: kg     
mwt_OW = 15.99940*1.6605E-27;  % Mol. wt. of OW 
mwt_HW = 1.008*1.6605E-27;     % Mol. Wt of HW 
m_water = mwt_water/NA;            % mass of 1 water molecule: kg     
m_OW = mwt_OW;                  % mass of OW atom: kg 
m_HW = mwt_HW;                  % mass of HW atom:kg 
beta = 1/(kb*T);                    % used in w_s_HO calculation
h = 6.626068E-34;                  % m^2 kg s^-1
N_water = 128;
int_md = 0.004;

mass(1) = m_OW;
mass(2) = m_HW;
mass(3) = m_HW;

% Input Section
disp('Reading velocities')
v_x = textread('vel-x-0-100ps.txt'); % Contains v_x of all molecules
v_y = textread('vel-y-0-100ps.txt'); % Contains v_y of all molecules
v_z = textread('vel-z-0-100ps.txt'); % Contains v_z of all molecules
disp('Reading Coords')
r_x = textread('coord-x-0-100ps.txt'); % Contains r_x of all molecules
r_y = textread('coord-y-0-100ps.txt'); % Contains r_y of all molecules
r_z = textread('coord-z-0-100ps.txt'); % Contains r_z of all molecules

v_x = v_x(:,1:385);
v_y = v_y(:,1:385);
v_z = v_z(:,1:385);

r_x = r_x(:,1:385);
r_y = r_y(:,1:385);
r_z = r_z(:,1:385);



time = v_x(:,1);

[m,n] = size (time);

disp('COM calculation')

%   N_water = 1; m = 1000;

  r_new_rx = zeros(m,3*N_water+1);
  r_new_x = zeros(m,3*N_water+1);
  
for g = 0:N_water-1
    
    for h = 1:m % for each time point
        
        I = zeros(3);
        r_comold_x(h,g+1) = 0;
        r_comold_y(h,g+1) = 0;
        r_comold_z(h,g+1) = 0;
        
        for l = 2:4 % for each molecule (3 atoms)


     % COM Calculation

            r_comold_x(h,g+1) = r_comold_x(h,g+1)+(mass(l-1)*(r_x(h,3*g+l)*(1E-9)));
            r_comold_y(h,g+1) = r_comold_y(h,g+1)+(mass(l-1)*(r_y(h,3*g+l)*(1E-9)));
            r_comold_z(h,g+1) = r_comold_z(h,g+1)+(mass(l-1)*(r_z(h,3*g+l)*(1E-9)));

        
        end
        
               
       r_comold_x(h,g+1) = r_comold_x(h,g+1)/(mass(1)+mass(2)+mass(3));
       r_comold_y(h,g+1) = r_comold_y(h,g+1)/(mass(1)+mass(2)+mass(3));
       r_comold_z(h,g+1) = r_comold_z(h,g+1)/(mass(1)+mass(2)+mass(3));
  
    end 
end

disp('Translating with respect to COM')

for re = 0:N_water-1
    for rf = 1:m
        for rg = 2:4

             r_new_rx(rf,3*re+rg) = r_comold_x(rf,re+1)-(r_x(rf,3*re+rg)*(1E-9));
             r_new_ry(rf,3*re+rg) = r_comold_y(rf,re+1)-(r_y(rf,3*re+rg)*(1E-9));
             r_new_rz(rf,3*re+rg) = r_comold_z(rf,re+1)-(r_z(rf,3*re+rg)*(1E-9));
        end
    end
end

disp( 'Calculating the moment of inertia tensor and principal axes of reference')

for le = 0:N_water-1
    for lf = 1:m
        
        I = zeros(3);
        
         for lg = 2:4

           I(1,1) = I(1,1) + mass(lg-1)*((r_new_ry(lf,3*le+lg))^2+(r_new_rz(lf,3*le+lg))^2);
           I(1,2) = I(1,2) - mass(lg-1)*(r_new_rx(lf,3*le+lg))*(r_new_ry(lf,3*le+lg));
           I(1,3) = I(1,3) - mass(lg-1)*(r_new_rx(lf,3*le+lg))*(r_new_rz(lf,3*le+lg));

           I(2,1) = I(2,1) - mass(lg-1)*(r_new_ry(lf,3*le+lg))*(r_new_rx(lf,3*le+lg));
           I(2,2) = I(2,2) + mass(lg-1)*((r_new_rx(lf,3*le+lg))^2+(r_new_rz(lf,3*le+lg))^2);
           I(2,3) = I(2,3) - mass(lg-1)*(r_new_ry(lf,3*le+lg))*(r_new_rz(lf,3*le+lg));

           I(3,1) = I(3,1) - mass(lg-1)*(r_new_rz(lf,3*le+lg))*(r_new_rx(lf,3*le+lg));
           I(3,2) = I(3,2) - mass(lg-1)*(r_new_rz(lf,3*le+lg))*(r_new_ry(lf,3*le+lg));
           I(3,3) = I(3,3) + mass(lg-1)*((r_new_rx(lf,3*le+lg))^2 + (r_new_ry(lf,3*le+lg))^2);
           
         end

        
        [V,D] = eig(I);
        
        D_vec = [D(1,1) D(2,2) D(3,3)];
        D_sort = sort(D_vec);
        
        if D_sort(1) == D(1,1)
            p1 = V(:,1);
        elseif D_sort(1) == D(2,2)
            p1 = V(:,2);
        else 
            p1 = V(:,3);
        end
        
        
        if D_sort(2) == D(1,1)
            p2 = V(:,1);
        elseif D_sort(2) == D(2,2)
            p2 = V(:,2);
        else 
            p2 = V(:,3);
        end
        
        
        if D_sort(3) == D(1,1)
            p3 = V(:,1);
        elseif D_sort(3) == D(2,2)
            p3 = V(:,2);
        else 
            p3 = V(:,3);
        end
        
        D1(lf,le+1) = D_sort(1);
        D2(lf,le+1) = D_sort(2);
        D3(lf,le+1) = D_sort(3);
                     
        p1_x(lf,le+1) = p1(1);
        p1_y(lf,le+1) = p1(2);
        p1_z(lf,le+1) = p1(3);
                
        p2_x(lf,le+1) = p2(1);
        p2_y(lf,le+1) = p2(2);
        p2_z(lf,le+1) = p2(3);
        
        p3_x(lf,le+1) = p3(1);
        p3_y(lf,le+1) = p3(2);
        p3_z(lf,le+1) = p3(3);
        
    end
end

% Angle between O-COM, H12 vectors and p1, p2, p3 axes 

% for he = 0:N_water-1
%     for hf = 1:m
%                     
%         p1_vec = [p1_x(hf,he+1) p1_y(hf,he+1) p1_z(hf,me+1)]';
%         p2_vec = [p2_x(hf,me+1) p2_y(hf,me+1) p2_z(hf,me+1)]';
%         p3_vec = [p3_x(hf,me+1) p3_y(hf,me+1) p3_z(hf,me+1)]';
%         
%         
%         
%         O_COM_old = [r_x(hf,3*he+2)*(1E-9)-r_comold_x(hf,he+1) r_y(hf,3*he+2)*(1E-9)-r_comold_y(hf,he+1) r_z(hf,3*he+2)*(1E-9)-r_comold_z(hf,he+1)]';
%         H12_old = [r_x(hf,3*he+3)-r_x(hf,3*he+4) r_y(hf,3*he+3)-r_y(hf,3*he+4) r_z(hf,3*he+3)-r_z(hf,3*he+4)]';
% 
%         p1_O_comold(hf,he+1) = acosd(dot(p1_vec,O_COM_old)/(norm(p1_vec)*norm(O_COM_old)));
%         p2_O_comold(hf,he+1) = acosd(dot(p2_vec,O_COM_old)/(norm(p2_vec)*norm(O_COM_old)));
%         p3_O_comold(hf,he+1) = acosd(dot(p3_vec,O_COM_old)/(norm(p3_vec)*norm(O_COM_old)));
% 
%         p1_H12old(hf,he+1) = acosd(dot(p1_vec,H12_old)/(norm(p1_vec)*norm(H12_old)));
%         p2_H12old(hf,he+1) = acosd(dot(p2_vec,H12_old)/(norm(p2_vec)*norm(H12_old)));
%         p3_H12old(hf,he+1) = acosd(dot(p3_vec,H12_old)/(norm(p3_vec)*norm(H12_old)));
% 
%     end 
% end



disp('Rotation Matrix Calculation and change of velocity and coordinate vectors')

x_axes = [1 0 0]';
y_axes = [0 1 0]';
z_axes = [0 0 1]'; 
   
    for vk = 0:N_water-1
        
        for gk = 1:m

        p1_vec = [p1_x(gk,vk+1) p1_y(gk,vk+1) p1_z(gk,vk+1)]';
        p2_vec = [p2_x(gk,vk+1) p2_y(gk,vk+1) p2_z(gk,vk+1)]';
        p3_vec = [p3_x(gk,vk+1) p3_y(gk,vk+1) p3_z(gk,vk+1)]';
        
        R(1,1) = dot(x_axes,p1_vec)/(norm(x_axes)*norm(p1_vec));
        R(1,2) = dot(y_axes,p1_vec)/(norm(y_axes)*norm(p1_vec));
        R(1,3) = dot(z_axes,p1_vec)/(norm(z_axes)*norm(p1_vec));
        
        R(2,1) = dot(x_axes,p2_vec)/(norm(x_axes)*norm(p2_vec));
        R(2,2) = dot(y_axes,p2_vec)/(norm(y_axes)*norm(p2_vec));
        R(2,3) = dot(z_axes,p2_vec)/(norm(z_axes)*norm(p2_vec));
        R(3,1) = dot(x_axes,p3_vec)/(norm(x_axes)*norm(p3_vec));
        R(3,2) = dot(y_axes,p3_vec)/(norm(y_axes)*norm(p3_vec));
        R(3,3) = dot(z_axes,p3_vec)/(norm(z_axes)*norm(p3_vec));
        
        detR(gk,vk+1) = det(R);
        
        if (det(R)) < 0 
         
            Rk = R*(-1);
            R = Rk;
           detR(gk,vk+1) = det(Rk);
        end
        
        
       for e = 2:4
           

       
           % Rotate coordinates to new axes
             
             r_new_x(gk,3*vk+e) = r_new_rx(gk,3*vk+e)*R(1,1) + r_new_ry(gk,3*vk+e)*R(1,2) + r_new_rz(gk,3*vk+e)*R(1,3);
             r_new_y(gk,3*vk+e) = r_new_rx(gk,3*vk+e)*R(2,1) + r_new_ry(gk,3*vk+e)*R(2,2) + r_new_rz(gk,3*vk+e)*R(2,3);
             r_new_z(gk,3*vk+e) = r_new_rx(gk,3*vk+e)*R(3,1) + r_new_ry(gk,3*vk+e)*R(3,2) + r_new_rz(gk,3*vk+e)*R(3,3); 

%              
             % Velocity vector transformation
             v_new_x(gk,3*vk+e) = v_x(gk,3*vk+e)*1E3*R(1,1) + v_y(gk,3*vk+e)*1E3*R(1,2) + v_z(gk,3*vk+e)*1E3*R(1,3);
             v_new_y(gk,3*vk+e) = v_x(gk,3*vk+e)*1E3*R(2,1) + v_y(gk,3*vk+e)*1E3*R(2,2) + v_z(gk,3*vk+e)*1E3*R(2,3);
             v_new_z(gk,3*vk+e) = v_x(gk,3*vk+e)*1E3*R(3,1) + v_y(gk,3*vk+e)*1E3*R(3,2) + v_z(gk,3*vk+e)*1E3*R(3,3);
  
             % Rotate principal axes
             
             p1_new_x(gk,vk+1) = p1_x(gk,vk+1)*R(1,1) + p1_y(gk,vk+1)*R(1,2) + p1_z(gk,vk+1)*R(1,3);
             p1_new_y(gk,vk+1) = p1_x(gk,vk+1)*R(2,1) + p1_y(gk,vk+1)*R(2,2) + p1_z(gk,vk+1)*R(2,3);
             p1_new_z(gk,vk+1) = p1_x(gk,vk+1)*R(3,1) + p1_y(gk,vk+1)*R(3,2) + p1_z(gk,vk+1)*R(3,3);
             
             p2_new_x(gk,vk+1) = p2_x(gk,vk+1)*R(1,1) + p2_y(gk,vk+1)*R(1,2) + p2_z(gk,vk+1)*R(1,3);
             p2_new_y(gk,vk+1) = p2_x(gk,vk+1)*R(2,1) + p2_y(gk,vk+1)*R(2,2) + p2_z(gk,vk+1)*R(2,3);
             p2_new_z(gk,vk+1) = p2_x(gk,vk+1)*R(3,1) + p2_y(gk,vk+1)*R(3,2) + p2_z(gk,vk+1)*R(3,3);
             
             p3_new_x(gk,vk+1) = p3_x(gk,vk+1)*R(1,1) + p3_y(gk,vk+1)*R(1,2) + p3_z(gk,vk+1)*R(1,3);
             p3_new_y(gk,vk+1) = p3_x(gk,vk+1)*R(2,1) + p3_y(gk,vk+1)*R(2,2) + p3_z(gk,vk+1)*R(2,3);
             p3_new_z(gk,vk+1) = p3_x(gk,vk+1)*R(3,1) + p3_y(gk,vk+1)*R(3,2) + p3_z(gk,vk+1)*R(3,3);
             
             
             
       end
        
        end
    end
        
   
    
% COM Calculation

for a = 0: N_water-1
    
   for b = 1:m
   
       r_com_x(b,a+1) = 0;
       r_com_y(b,a+1) = 0;
       r_com_z(b,a+1) = 0;
       
       for c=2:4
            
           r_com_x(b,a+1) = r_com_x(b,a+1)+(mass(c-1)*(r_new_rx(b,3*a+c)));
           r_com_y(b,a+1) = r_com_y(b,a+1)+(mass(c-1)*(r_new_ry(b,3*a+c)));
           r_com_z(b,a+1) = r_com_z(b,a+1)+(mass(c-1)*(r_new_rz(b,3*a+c)));
           
       end
       
       r_com_x(b,a+1) = r_com_x(b,a+1)/(mass(1)+mass(2)+mass(3));
       r_com_y(b,a+1) = r_com_y(b,a+1)/(mass(1)+mass(2)+mass(3));
       r_com_z(b,a+1) = r_com_z(b,a+1)/(mass(1)+mass(2)+mass(3));
       
        
   end
          
end

% Angle between O-COM, H12 vectors and p1, p2, p3 axes 

for me = 0:N_water-1
    for mf = 1:m
                    
        p1_vec = [p1_new_x(mf,me+1) p1_new_y(mf,me+1) p1_new_z(mf,me+1)]';
        p2_vec = [p2_new_x(mf,me+1) p2_new_y(mf,me+1) p2_new_z(mf,me+1)]';
        p3_vec = [p3_new_x(mf,me+1) p3_new_y(mf,me+1) p3_new_z(mf,me+1)]';
        
        O_COM = [r_new_x(mf,3*me+2)-r_com_x(mf,me+1) r_new_y(mf,3*me+2)-r_com_y(mf,me+1) r_new_z(mf,3*me+2)-r_com_z(mf,me+1)]';
        H12 = [r_new_x(mf,3*me+3)-r_new_x(mf,3*me+4) r_new_y(mf,3*me+3)-r_new_y(mf,3*me+4) r_new_z(mf,3*me+3)-r_new_z(mf,3*me+4)]';

        p1_O_com(mf,me+1) = acosd(dot(p1_vec,O_COM)/(norm(p1_vec)*norm(O_COM)));
        p2_O_com(mf,me+1) = acosd(dot(p2_vec,O_COM)/(norm(p2_vec)*norm(O_COM)));
        p3_O_com(mf,me+1) = acosd(dot(p3_vec,O_COM)/(norm(p3_vec)*norm(O_COM)));

        p1_H12(mf,me+1) = acosd(dot(p1_vec,H12)/(norm(p1_vec)*norm(H12)));
        p2_H12(mf,me+1) = acosd(dot(p2_vec,H12)/(norm(p2_vec)*norm(H12)));
        p3_H12(mf,me+1) = acosd(dot(p3_vec,H12)/(norm(p3_vec)*norm(H12)));

    end 
end
% 
% 
% % Distance between the atoms in the new coordinates
% 
% for fg = 0:N_water-1
%     for gf = 1:m
%                  
%             r_rnew_OH1(gf,fg+1) = sqrt((r_new_rx(gf,3*fg+2)-r_new_rx(gf,3*fg+3))^2 + (r_new_ry(gf,3*fg+2)-r_new_ry(gf,3*fg+3))^2 +(r_new_rz(gf,3*fg+2)-r_new_rz(gf,3*fg+3))^2);
%             r_rnew_OH2(gf,fg+1) = sqrt((r_new_rx(gf,3*fg+2)-r_new_rx(gf,3*fg+4))^2 + (r_new_ry(gf,3*fg+2)-r_new_ry(gf,3*fg+4))^2 +(r_new_rz(gf,3*fg+2)-r_new_rz(gf,3*fg+4))^2);
%             
%             r_new_OH1(gf,fg+1) = sqrt((r_new_x(gf,3*fg+2)-r_new_x(gf,3*fg+3))^2 + (r_new_y(gf,3*fg+2)-r_new_y(gf,3*fg+3))^2 +(r_new_z(gf,3*fg+2)-r_new_z(gf,3*fg+3))^2);
%             r_new_OH2(gf,fg+1) = sqrt((r_new_x(gf,3*fg+2)-r_new_x(gf,3*fg+4))^2 + (r_new_y(gf,3*fg+2)-r_new_y(gf,3*fg+4))^2 +(r_new_z(gf,3*fg+2)-r_new_z(gf,3*fg+4))^2);
%            
%             r_OH1(gf,fg+1) = sqrt((r_x(gf,3*fg+2)-r_x(gf,3*fg+3))^2 + (r_y(gf,3*fg+2)-r_y(gf,3*fg+3))^2 +(r_z(gf,3*fg+2)-r_z(gf,3*fg+3))^2);
%             r_OH2(gf,fg+1) = sqrt((r_x(gf,3*fg+2)-r_x(gf,3*fg+4))^2 + (r_y(gf,3*fg+2)-r_y(gf,3*fg+4))^2 +(r_z(gf,3*fg+2)-r_z(gf,3*fg+4))^2);
%             
%             r_new_H12(gf,fg+1) = sqrt((r_new_x(gf,3*fg+3)-r_new_x(gf,3*fg+4))^2 + (r_new_y(gf,3*fg+3)-r_new_y(gf,3*fg+4))^2 +(r_new_z(gf,3*fg+3)-r_new_z(gf,3*fg+4))^2);
%             
%             r_Ocom(gf,fg+1) = sqrt((r_new_x(gf,3*fg+2)-r_new_x(gf,fg+1))^2 + (r_new_y(gf,3*fg+2)-r_new_y(gf,fg+1))^2 +(r_new_z(gf,3*fg+2)-r_com_z(gf,fg+1))^2);
%             r_H1com(gf,fg+1) = sqrt((r_new_x(gf,3*fg+3)-r_new_x(gf,fg+1))^2 + (r_new_y(gf,3*fg+3)-r_new_y(gf,fg+1))^2 +(r_new_z(gf,3*fg+3)-r_com_z(gf,fg+1))^2);
%             r_H2com(gf,fg+1) = sqrt((r_new_x(gf,3*fg+4)-r_new_x(gf,fg+1))^2 + (r_new_y(gf,3*fg+4)-r_new_y(gf,fg+1))^2 +(r_new_z(gf,3*fg+4)-r_com_z(gf,fg+1))^2);
%          
%     end 
% end
         
disp('Angular momentum calculation')

for k = 0: N_water-1
   
    for j = 1:m
        
        L =0;
       
       for l = 2:4 % foreach molecule (3 atoms)
     
             r_vec = [(r_new_x(j,3*k+l)) r_new_y(j,3*k+l) r_new_z(j,3*k+l)]; 
        
             v_vec = [v_new_x(j,3*k+l) v_new_y(j,3*k+l) v_new_z(j,3*k+l)];
             
             mag_vnew(j,3*k+l) = sqrt(v_vec(1)^2+v_vec(2)^2+v_vec(3)^2);
             mag_v(j,3*k+l) = sqrt(v_x(j,3*k+l)^2+v_y(j,3*k+l)^2+v_x(j,3*k+l)^2);
             
             L = L + (mass(l-1) *(cross(r_vec,v_vec)));
       end
       
       L_x(j,k+1) = L(1);
       L_y(j,k+1) = L(2);
       L_z(j,k+1) = L(3);

% moment of inertia tensor
        
        I1 = zeros(3);
        
        for z = 2:4 % for each molecule (3 atoms)
%           
             
            
            I1(1,1) = I1(1,1) + mass(z-1)*((r_new_y(j,3*k+z))^2+(r_new_z(j,3*k+z))^2);
            I1(1,2) = I1(1,2) + mass(z-1)*(r_new_x(j,3*k+z))*(r_new_y(j,3*k+z));
            I1(1,3) = I1(1,3) + mass(z-1)*(r_new_x(j,3*k+z))*(r_new_z(j,3*k+z));

            I1(2,1) = I1(2,1) + mass(z-1)*(r_new_y(j,3*k+z))*(r_new_x(j,3*k+z));
            I1(2,2) = I1(2,2) + mass(z-1)*((r_new_x(j,3*k+z))^2+(r_new_z(j,3*k+z))^2);
            I1(2,3) = I1(2,3) + mass(z-1)*(r_new_y(j,3*k+z))*(r_new_z(j,3*k+z));

            I1(3,1) = I1(3,1) + mass(z-1)*(r_new_z(j,3*k+z))*(r_new_x(j,3*k+z));
            I1(3,2) = I1(3,2) + mass(z-1)*(r_new_z(j,3*k+z))*(r_new_y(j,3*k+z));
            I1(3,3) = I1(3,3) + mass(z-1)*((r_new_x(j,3*k+z))^2 + (r_new_y(j,3*k+z))^2);

                  
        end
        
       
         
        I_x(j,k+1) = I1(1,1);
        I_y(j,k+1) = I1(2,2);
        I_z(j,k+1) = I1(3,3);
        
        I_xy(j,k+1) = -I1(1,2);
        I_xz(j,k+1) = -I1(1,3);
                
        I_yx(j,k+1) = -I1(2,1);
        I_yz(j,k+1) = -I1(2,3);
        
        I_zx(j,k+1) = -I1(3,1);
        I_zy(j,k+1) = -I1(3,2);
                                     
        w_vec = inv(I1)*L'; 
        w(j,k+1) = norm(w_vec);
        w_x(j,k+1) = w_vec(1);
        w_y(j,k+1) = w_vec(2);
        w_z(j,k+1) = w_vec(3);
        
        ke_tot(j,k+1) = (I1(1,1)*(w_x(j,k+1)^2) + I1(2,2)*(w_y(j,k+1)^2) + I1(3,3)*(w_z(j,k+1)^2))*0.5;
        
        ke_x(j,k+1) = 0.5*(I1(1,1)*(w_x(j,k+1)^2));
        ke_y(j,k+1) = 0.5*(I1(2,2)*(w_y(j,k+1)^2));
        ke_z(j,k+1) = 0.5*(I1(3,3)*(w_z(j,k+1)^2));

    
    end       
end

disp(' Autocorrelation Calculations')

% Iy1 = zeros(1,12500); Iy2 = zeros(1,12500); Iy3 =zeros(1,12500);C_Rot = 0;

Iy1 = 0;Iy2 = 0;Iy3 = 0;C_Rot = 0;

% Block Averaging
% Divide w_x into blocks and calculate the ACF for each. 
% 
% blocksize = 10ps ; I have 100 ps data. 

nblock = m*int_md; % ps

tb = nblock*(1/int_md);
block_num = round(m/tb);

 for de = 0:N_water-1
     for bn = 0:block_num-1
         
         de
         bn
         
         meanI_x = mean(I_x(bn*tb+1:(bn+1)*tb,de+1));
         meanI_y = mean(I_y(bn*tb+1:(bn+1)*tb,de+1));
         meanI_z = mean(I_z(bn*tb+1:(bn+1)*tb,de+1));

         y1 = auto(w_x(bn*tb+1:(bn+1)*tb,de+1));
         y2 = auto(w_y(bn*tb+1:(bn+1)*tb,de+1));
         y3 = auto(w_z(bn*tb+1:(bn+1)*tb,de+1));

         Iy1 = Iy1+meanI_x*y1;
         Iy2 = Iy2+meanI_y*y2;
         Iy3 = Iy3+meanI_z*y3;

    end
 end

 %   
C_x = Iy1/(N_water*block_num);
C_y = Iy2/(N_water*block_num);
C_z = Iy3/(N_water*block_num);

C_Rot = C_x + C_y + C_z;% unit: (rad/s)^2
% C_Rot = C_Rot/(4*pi*pi); % unit: (turns/s)^2

crot = C_Rot';


y1 = 0;y2 = 0;y3 = 0;Ny1 = 0;Ny2 = 0;Ny3 = 0;
% Non-weighted autocorrelations
 for ve = 0:N_water-1
     for vn = 0:block_num-1
                  
         y1 = auto(w_x(vn*tb+1:(vn+1)*tb,ve+1));
         y2 = auto(w_y(vn*tb+1:(vn+1)*tb,ve+1));
         y3 = auto(w_z(vn*tb+1:(vn+1)*tb,ve+1));

         Ny1 = Ny1+y1;
         Ny2 = Ny2+y2;
         Ny3 = Ny3+y3;
     end
 end


%   
NC_x = Ny1/(N_water*block_num);
NC_y = Ny2/(N_water*block_num);
NC_z = Ny3/(N_water*block_num);

NC_Rot = NC_x + NC_y + NC_z;

save ('ncrot.txt','NC_Rot','-ascii')
save ('c_rot.txt','crot','-ascii')
save ('wy_major.xvg','w_y','-ascii')
save ('wx_major.xvg','w_x','-ascii')
save ('wz_major.xvg','w_z','-ascii')

save ('wy_major.txt','w_y','-ascii')
save ('wx_major.txt','w_x','-ascii')
save ('wz_major.txt','w_z','-ascii')

save('I_x_major.txt','I_x','-ascii')
save('I_y_major.txt','I_y','-ascii')
save('I_z_major.txt','I_z','-ascii')


% % Fourier Spectra
% % 
% t = time(1:(m+1)/2);
% [gt0, area] = symmfft(t,C_Rot)
% %         

% Rotational Tempertaure Calculations

% B = h/(8 π c I)
% 
% Where B is the rotational constant (cm-1) 
% h is Plancks constant (gm cm2/sec) 
% c is the speed of light (cm/sec) 
% I is the moment of inertia (gm cm2)

h_cm = 6.626068E-27;  % g cm^2/s
c_cm = 3*1E10; % cm/s
I_convfac = 1E7; % kg m^2 --> g cm^2

for ie = 0:N_water-1
    for je = 1:m
        
        B_cm_x(je,ie+1) = (h_cm)/(8*pi*pi*c_cm*I_x(je,ie+1)*I_convfac);
        B_cm_y(je,ie+1) = (h_cm)/(8*pi*pi*c_cm*I_y(je,ie+1)*I_convfac);
        B_cm_z(je,ie+1) = (h_cm)/(8*pi*pi*c_cm*I_z(je,ie+1)*I_convfac);
        
        theta_x(je,ie+1) = B_cm_x(je,ie+1)/0.6950;  % kb = 0.695 cm^-1 K^-1
        theta_y(je,ie+1) = B_cm_y(je,ie+1)/0.6950;
        theta_z(je,ie+1) = B_cm_z(je,ie+1)/0.6950;
        
    end
end


save ('wy_major.xvg','w_y','-ascii')
save ('wx_major.xvg','w_x','-ascii')
save ('wz_major.xvg','w_z','-ascii')

save ('wy_major.txt','w_y','-ascii')
save ('wx_major.txt','w_x','-ascii')
save ('wz_major.txt','w_z','-ascii')

save('I_x_major.txt','I_x','-ascii')
save('I_y_major.txt','I_y','-ascii')
save('I_z_major.txt','I_z','-ascii')

